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Abstract. We develop a variational scheme in a field theoretic approach to a 
stochastic process. While various stochastic processes can be expressed using master 
equations, in general it is difficult to solve the master equations exactly, and it is also 
hard to solve the master equations numerically because of the curse of dimensionality. 
The field theoretic approach has been used in order to study such complicated 
master equations, and the variational scheme achieves tremendous reduction in the 
dimensionality of master equations. For the variational method, only the Poisson 
ansatz has been used, in which one restricts the variational function to a Poisson 
distribution. Hence, one has dealt with only restricted fluctuation effects. We develop 
the variational method further, which enables us to treat an arbitrary variational 
function. It is shown that the variational scheme developed gives a quantitatively 
good approximation for master equations which describe a stochastic gene regulatory 
network. 
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1. Introduction 

Master equations describe various stochastic phenomena. For example, a reaction- 
diffusion process, which is one of the examples of non-equilibrium systems, is expressed 
by a master equation. Usually, it is difficult to obtain the exact solution of the master 
equation because of the non-linearity of the corresponding moment equations or its high 
dimensionality. The direct numerical solution is also difficult to obtain because there are 
an enormous number of coupled equations to be solved. While numerical simulations 
such as the Gillespie algorithm [1] are available for studying complicated stochastic 
systems, a coarse-grained analytical approach would be more worthwhile. The field 
theoretic approach to the reaction-diffusion process has achieved significant successes [2] . 
The analogy of the master equation to a quantum system has been introduced by 
Doi [3,4], and several authors revived the formalism [5,6]. The field theoretic approach 
has revealed the anomalous kinetics in the react ion- diffusion systems incorporating the 
renormalization group method [2]. In addition, the field theoretic description has been 
used not only for the reaction-diffusion processes, but also for various phenomena such 
as packet flow [7], the Malthus-Verhulst process [8], stochastic sandpile models [9,10], 
and neural networks [11]. 

Recently, Sasai and Wolynes [12] have developed the field theoretic approach to 
a stochastic gene network. The gene network consists of active and inactive genes, 
proteins produced by the genes, and a mechanism of switching between the active and 
inactive states caused by the regulatory proteins. The complicated system is described 
by a set of master equations, as in the case with the reaction-diffusion process. For the 
case of only one gene, the exact solution has been obtained [13], but when one consider a 
general case, i.e., a gene regulatory network, it is difficult to solve the master equations 
exactly. We therefore need an approximation method. Sasai and Wolynes [12] have 
used the variational method for non-equilibrium systems which has been proposed by 
Eyink [14,15], The variational method gives us an efficient scheme of approximation for 
complicated master equations; it can achieve enormous reduction in the dimensionality 
of the problem by solving variationally the quantum field theoretic equations which 
are obtained by the original master equations. This means that the variational scheme 
reduces the coupled master equations with a huge number of variables to a set of ordinary 
differential equations with a small number of parameters. 

So far, several schemes of approximation for master equations have been proposed 
[16, 17]. In the system size expansion or the Kramers- Moyal expansion [17], master 
equation with 'discrete' variables are substituted into a Fokker-Planck equation with 
'continuous' variables. While the differential equation with continuous variables is easier 
to treat, these approximation schemes are applicable only for a system with the large 
size. The moment equation approach [17] can be used even for small systems, and 
it gives an exact solution when the moment equations are closed. However, if the 
moment equations have non-linear terms and are not closed, it is difficult to obtain 
the exact solution. To our knowledge, there is no systematic scheme which produces a 
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closed set of moment equations. The variational scheme gives a closed set of equations 
in a systematic way, and therefore the variational scheme is expected as a candidate 
systematic approximation scheme for complicated master equations. The variational 
scheme in [12,18,19] is based on the Poisson ansatz, in which the mean and the variance 
of the variational function are the same. It has been revealed that the solutions obtained 
by using the Poisson ansatz are correct only qualitatively for the repressilator system 
with two genes [18,19]. 

The aim of the present paper is to develop a variational scheme beyond the 
Poisson ansatz. In principle, the variational function should be a discrete probability 
distribution. The Poisson distribution (the Poisson ansatz) is useful for the functional 
variation because the Poisson ansatz corresponds to the coherent state in the field 
theoretic description. On the other hands, the other discrete probability distribution 
is difficult to treat in the variational scheme. In order to avoid the difficulty of the 
variational calculations, we propose the use of the superposition of the coherent states 
as the variational function. By using the superposition of the coherent states, it becomes 
possible to assume an arbitrary continuous probability function as the variational 
function. We will apply the variational method to a gene regulatory network, which 
is the same as the one in [18], and confirm that the new method gives a quantitatively 
correct solution. 

The construction of the present paper is as follows. In section 2, we define the 
gene regulatory network and master equations to be solved, and we also give the field 
theoretic description of the master equations. The variational scheme proposed by 
Sasai and Wolynes are reviewed in section 3. Section 4 is the main part of the present 
paper, and gives the new variational function beyond the Poisson ansatz. The numerical 
experiments are also performed in order to confirm the validity of the new scheme. 
Finally, we give some concluding remarks in section 5. 

2. Model and formalism 

2.1. Master equations of a gene regulatory network 

We here give an explicit example, i.e., a chemical reaction network involved in gene 
regulations, which has been used in [12] and [18]. The master equations for the gene 
regulatory network give a closed set of moment equations, and hence we can confirm 
the validity of the variational scheme by comparing the results of the variational scheme 
with those obtained by the moment equations. 

Figure [1] shows the gene regulatory network. In the network, there are two genes 
which are labeled by A and B, respectively. Each gene produces a repressor protein 
which binds to the operator site of the other gene changing the activity. When gene 
a (a = A or B) is not bound by the repressor proteins, the gene can produce its own 
proteins at the rate g a \. The gene bound by the repressor proteins produces own proteins 
at the rate g a0 . Note that the subscripts '1' and '0' represent the active state and the 
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Figure 1. Illustration of the gene regulatory network. The Protein produced by gene 
A is a repressor which binds to gene B, and vice versa. The production rate g depends 
on whether the gene is bound or not bound by the repressor. 



inactive state of the gene, respectively. Each protein spontaneously degrades, and the 
degradation rate is k a . The rate of binding of the proteins to a gene and the rate of 
detaching from a gene are represented by h a and f a , respectively. In the present paper, 
we consider the case where dimer proteins repress the expression of a gene. 

The next step is to write down master equations for the gene regulatory network. 
Hereafter, the number of proteins produced by gene a is denoted as n a . Using the two 
component vector notation 

3 (£&:!))• (1) 

the master equation for the probability with which there are n a proteins can be written 

as 

^P a {n w t)=( 9 f ° )[P a (n a -l,t)-P a (n on t)] 
at \ U g a0 ) 

+k a {(n a + l)P a (n Q + 1, t) - n a P a (n a , t)) 

where (a,/?) = {(A, B), (B, A)}. Although we might be able to solve the 
master equation ([2]) numerically, it becomes difficult to solve the master equation 
numerically when the number of genes increases. Even for only one gene, we have 
2 x (the number of state n) coupled differential equations. In order to reduce the 
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dimensionality of the problem, we use the field theoretic description and a variational 
scheme. 

2.2. Field theoretic description 

It is revealed that the quantum field theoretic method is useful to solve the master 
equations. We briefly review the quantum field theoretic description for the gene 
regulatory network [12]. 

First of all, we define the ket vector \n) as the state in which there are n proteins in 
the system. For each protein (protein A and protein B), a creation and an annihilation 
operators are introduced: 



al\n a ) 



n a + 1), 



(3) 



a a \n a ) = n a \n a - 1), (4) 

where the index a is A or B. The creation and the annihilation operators satisfy the 
following commutation relation 

[a a ,4] = l, (5) 

and the vacuum state |0 a ) and its conjugate (0 a | are defined to satisfy 




Note that the n-proteins state |n) is not normalized in the usual sense, but the states 
are orthogonal, because (n\m) = m\5 n ^ mi where 5 n ^ m is the Kronecker delta, 

Using the above quantum field theoretic formalism, we write the state which 
corresponds to a probability distribution vector P a (n a ,t) as 




(8) 



The state \ip a ) only describes the state of gene a, and hence the state of the whole 
system is denoted by 






When the total Hamiltonian operator Vt is defined as 
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the Hamiltonian fi a which operates only gene a is derived from the original master 
equation as 

Q = ( g a i(al - 1) + k a (a a - aia a ) ^ ® 1 

( fa \ , , ( =M4) 2 (a*) 2 \ 

where the suffix a or j3 of each operator means that the operator acts only on gene a 
or (3. The first term corresponds to the birth-death part of proteins, and plays a role in 
the diffusion effects. The second and third terms represent the interactions between two 
genes. Note that the "Hamiltonian" is non-Hermitian and it is a little different from 
the ordinary quantum mechanics one. For instances, expected values are linear not 
bilinear in \ip a ), and averages for \ip a ) are obtained by taking the scalar product with 
the bra vector ((0|e a " (0|e a "). However, in spite of the non-Hermitian property and a 
slight difference from the ordinary quantum mechanics, many quantum field theoretic 
techniques may be applied, albeit with some modifications. 



3. Variational approach 

3.1. Variation of the effective action 

In order to reduce the dimensionality of the problem, a variational method developed 
by Eyink [14, 15] can be used. We here briefly review the method [12]. 
When we define an effective action V as 



dt($\(d t -n)\v), (13) 

equation ([TO]) is equivalent to the functional variation = 0. Because of the non- 

Hermitian property, it is not always true that the left eigenvectors and right eigenvectors 
are the same. Hence, we assume two variational functions for the bra and ket states, 
respectively. We assume that the ket state (or the bra state ($|) is parametrized 
by (or x L ), and where x fi and x L are vectors with K components: 

x iJ = {xf,x^,---,x|}, (14) 

X = { X l ? X 2 ? ' ' ' ' X I<}- (15) 

A set of finite dimensional equations for parameters x fi and x L is obtained by the 
functional variation procedure. Note that we set $(x L = 0) to be consistent with the 
probabilistic interpretation, so that 

($(x L = 0)|vl/(x fi )) = l, (16) 

which is the normalization condition for the probability distribution. We, therefore, 
obtain the following equation which stems from an extremum of the action: 



K 



y.(— 

^ \ dx L 

,1=1 \ u - h m 



dx? / dt \ dxk, 



0, for m= 1,2, (17) 
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The only remaining procedure is to give an explicit ansatz for ($| and 
This corresponds to the fact that we restrict a probability distribution P a i(n a ,t) (or 
Pao( n a,t)) to a specific form with a few free parameters. Although it is difficult to 
calculate the time evolution of the probability distribution directly, the variational 
scheme enables us to get a set of time evolution equations for the time-dependent 
parameters x^; the time evolution equations for the parameters are determined 
variationally through equation ffTTl) . Note that in the variational scheme, it is necessary 
to set an adequate ansatz for ($| and in order to get qualitatively or quantitatively 
correct results. 



3.2. Poisson ansatz 

As for the choice of the ansatz in equation (fTTl) . only the Poisson ansatz has been 
proposed so far [12, 18]. The Poisson ansatz is a reasonable choice because the 
steady-state probability distribution for a simple birth-death problem is the Poisson 
distribution. Furthermore, the Poisson ansatz is based on the coherent state, which 
makes it easy to perform the variational calculation. 

In the Poisson ansatz, we assume the following ket vector 




(18) 



C a i exp X al (a) a - 1 
^ C a0 exp [X a0 (al - 1 

and as the bra ansatz, 

(<f> a \ = ( (0«| exp (a a + Xai + A^iV) (0 a | exp (a a + a£q + A ao a a) ) • (19) 

Note that although one might construct a bra ansatz with a lot of free parameters, the 
same number of free parameters for the bra ansatz as that of the ket ansatz is enough 
for constructing the time evolution equations for the free parameters in the ket ansatz. 
We therefore have in total 16 parameters in the bra and ket variational functions; 

x^ = {Cai, Cao, Xai, Xao, Cbi, Cbo, X B i, X B q}, (20) 
Y i_r\(o) x (o) x (i) x (i) x (o) x (o) x (i) x (ih / 91 x 

Performing the variational calculation of equation (|T7|) . we finally have six coupled 
ordinary differential equations [18]; the number of parameters for the ket ansatz is eight 
but there are two constraints from the normalization of the probability: Cai + Cao — 1 
and Cbi + Cbo — 1- I n addition, all parameters in the bra ansatz are set to be zero 
finally, and therefore there are only six equations. 



4. Beyond the Poisson ansatz 

Although it has been shown that the Poisson ansatz gives qualitatively appropriate 
results for the gene regulatory network [18, 19], the solution of the Poisson ansatz is not 
quantitatively correct. Hence, it is necessary to develop the variational scheme beyond 
the Poisson ansatz. 
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In general, a state in the field theoretic description is described by J2^=o P( n )\ n )i 
where P(n) is a discrete probability distribution. Note that P(ri) must be a discrete 
probability distribution because n takes an integer value. When we use the Poisson 
distribution as the probability P(n), we have the coherent states and then it is easy 
to calculate the functional variation. However, for the other discrete probability 
distribution, it is difficult to calculate the functional variation in equation ffTTj) . 

In order to overcome the problems, we here propose a new ansatz for the variational 
scheme. The new ansatz is based on the idea in which we use the superposition of the 
coherent states. For example, when we want to have two parameters for the variational 
function, the following ansatz for the ket state should be used: 

u\ = ( C al J^dxF(x\^j,fiap^[x(al - l)]|0 a ) j ^ 
\ C a oJ^°dxF(x;fi ( ^,fi { ^)exp[x(al - l)]|0 a ) J 

The new ansatz, the superposition ansatz, means that we take a superposition of the 
Poisson distributions with different mean values. The 'continuous' variational function 
F(x) is a probability density with two parameters. In the gene regulatory networks, 
the state \n) does not have negative n, so the integral range of F(x) should be taken as 
x > 0. We note that the formalism can be extended to the case with more complicated 
variational function with many parameters. The ansatz with only two free parameters 
in equation fl22l) is a simple case beyond the Poisson ansatz. 

Using the superposition ansatz, we can easily perform the variational calculation 
because the variational functions are based on the coherent states. In addition, the 
superposition ansatz enables us to use a continuous variational function. Unlike using 
continuous approximations of master equations, such as the Kramers-Moyal expansion 
and the system size expansion [16,17], the use of the continuous variational function 
in the superposition ansatz includes the discrete characteristics of the original master 
equation due to the use of the coherent states. 

As the bra ansatz, we here simply take 

(0 a | exp (a a + AS + A$a a + Xai(a a )^ 



(0 a | exp [a a + A„i + Xao a a + A c (2) ' 




(23) 



x a0 1 

where T represents the transposed matrix. Finally, we have the following 24 parameters 
for the variational calculation 

R — fn n ( 1 ) W ( 2 ) ( 2 ) n n W W ( 2 ) ( 2 h (oa\ 

X — {L^Al, <^A0, /^AlJ A 1 AO' AfAl) ^AO? U B1, <^B0, ^BIi ^BQi ^B\i VboTj l 24 J 
v i_J\(0) \(0) x(2) x (2) ,(0) ,(0) ,(1) x (l) x (2) ,(2), ft)t ^ 

Using the superposition ansatz of equation (|22|) . we have 10 ordinary differential 
equations to be solved by numerical integration. (The ket ansatz has 12 parameters, 
but there are two constraints related to the normalization of the probability, so that we 
have only 10 equations.) 
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Figure 2. Probabilities Cax and Cbx in the long time limit (in the steady state). 
The horizontal axis means the rescaled parameter X &( ^ = (gi + go) /(2k a). At a certain 
critical point, there is the bifurcation from the monostable state to the bistable state. 
The values of Cax and Cbi are represented by the same symbol for simplicity. Wc 
note that Cai and Cbi take different stable states to each other in the bistable state. 
The initial state determines which state (Cai or Cbi) takes the higher value in the 
bistable state. 



In what follows, we check the superposition ansatz by numerical experiments. As 
the variational function with two parameters, we here take a gamma distribution; 

FfeM)=^'^y^. (26) 

The gamma function has the mean k9 and the variance k6 2 . As in the case of the 
Poisson ansatz, a set of ordinary differential equations for the parameters are obtained 
by using a simple symbolic algebraic calculation in the field theoretic description. The 
resulting equations are a little long, so we write the resulting equations in the appendix. 

We performed a numerical experiment in order to confirm the improvement achieved 
by the superposition ansatz. We fixed all parameters except the protein synthesis 
rate g x = g A1 = g m ; k A = k B = 1, f A = f B = 0.5, h A = h B = /a/500, and 
g = g A0 = g B o = 0, which are the same parameters as in [18]. For various initial 
states of the variational parameters, the steady state is obtained in the long time limit. 
Figure [2] shows the probabilities C A \ and C B \ with which genes A and B are in the 
active state, as a function of X ac j = (g% + g Q )/(2k A ). The values of C A \ and C B \ are 
represented by the same symbol for simplicity. As shown in figure [2], the bifurcation 
from the monostable state to the bistable state is observed as increasing X & ^. We note 
that in the monostable state the values of G A \ and C B \ are the same, but C A \ and C B \ 
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take different stable states to each other in the bistable state. It depends on the initial 
parameters which probability, Cai or Cbi, is larger than the other in the bistable state. 

The solid line in figure [2] is obtained from the moment equations in [18], which is a 
closed set of equations and gives exact solutions for the present case. The Poisson ansatz 
gives a qualitatively good results; the bifurcation is observed. However, the bifurcation 
point is different from the result from the moment equations. In contrast, the results 
from the gamma distribution ansatz are in quantitatively good agreement with the 
moment equations. The numerical results confirm the validity of the superposition 
ansatz. 

In the moment equations in reference [18], the first and second moments of the 
protein number need to be taken into account. In the Poisson ansatz, the mean and the 
variance should be same, and then the Poisson ansatz does not give the quantitatively 
correct results because the second moment depends on the first moment. In contrast, 
the gamma distribution ansatz includes two free parameters so the second moment 
of x is independent of the first moment. We consider that this is why the gamma 
distribution ansatz gives quantitatively correct results. In addition, we have checked 
that the other ansatz, e.g., a log normal distribution, also works well for calculating the 
bifurcation point correctly. However, we note that there may be a suitable variational 
function for investigating higher correlations in the gene regulatory networks. Study of 
the applicability of the variational scheme will be important in the future. 

5. Concluding remarks 

In the present paper, a new ansatz for the variational scheme was proposed. The 
superposition ansatz is based on the coherent states, so it gives us a straightforward 
extension of the variational scheme with the Poisson ansatz. In addition, it enables us to 
use various continuous probability densities as the variational function. The availability 
of the superposition ansatz was confirmed in a simple gene regulatory network. The 
superposition ansatz gives a quantitatively correct solution, while the Poisson ansatz is 
adequate only qualitatively. 

The concept of the superposition of the Poisson distributions seems to be related to 
the Poisson representation [17]. The coefficient function in the Poisson representation 
can take complex numbers, so that it is not always true that the coefficient function 
corresponds to the probability distribution. The relationship between the Poisson 
representation and the quantum field theoretic representation has been pointed out [20], 
and actually, our variational scheme is related to the Poisson representation; it is easy 
to see that the superposition ansatz restricts the coefficient function in the Poisson 
representation to being a certain variational function. This correspondence between the 
superposition ansatz and the Poisson representation would give us further extensions of 
the superposition ansatz; it might be possible to use a function of complex variable as 
the variational function. This is a future work. 

The variational method and the quantum field theoretical description would give 
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new and useful schemes of approximation for complicated master equations. For 
example, the superposition ansatz enables us to extend the variational scheme to 
multivariate cases [21]. These approximation methods are important for researching 
complex systems such as biological systems and social systems. Furthermore, it may 
be possible to study the complex systems more analytically by using the quantum 
field theoretical description. Such researches would give deep insight into the complex 
systems. 



Appendix A. Time evolution equations in the superposition Ansatz 



From equation (TTTI) and the superposition ansatz of (1221 . a set of coupled ordinary 
differential equations are derived. Here, we use the following notation for simplicity: 
F(x; Hal, A*ai ) = F a i(x). Performing the variational calculation, we obtain the following 
five time-evolution equations for the parameters related to gene A: 

dC Al 



dt 

dC A1 
dt Jo 



-C.i C 



(h a r°° h a r°° \ 

Cm YJo fa2F a( s ) + ^y/ o dxx 2 F B0 (x)\+f A C A0 , (A.r 



dC A0 
dt 



dC A1 
dt Jo 



dxxF A1 (x) + C A1 



C A \ 



dt Jo 

fOO 



dt Jo 



°°, dF A1 (x) 
axx- 



9A1 



poo 

+ C A0 f A / dxxF A0 (x) 
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dxF A i (x 
dx A x A F A1 (x A ) 
x \C m dx B x 2 B F BX {x B ) + C B0 J dx B x 2 B F B0 (x B 

dF A0 (x 



dxxFA n (x) + C 



d^ A l 

,_ dT 



dxx- 



d/j, 



(i) 

.40 



+ C 



A ( 2 ) 



°°, dF A0 (x) 
axx- 



djj. 



(2) 
.40 



c 



AO 



fOO poo 

g A0 -k dxF A0 (x) -C A0 f A dxxF A0 (x) 
Jo J Jo 



h A f°° 
+— Cm / dx A x A F A1 (x A ) 
2 Jo 



x\C B1 dx B x B F B1 (x B ) + C B0 / dx B x B F B0 (x B ) 



dxx 2 F A1 (x) + C A1 d ^ M 



« dF A1 (x' A - {2) 



(A.2) 



(A.3) 



dt Jo 



dxx z 



C A \ 
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axx 
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dxF Al (x) 2g A1 x-2kx 2 + C A0 f A dxx 2 F A0 (x 
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(A.4) 



A field theoretic approach to master equations and a variational method beyond the Poisson ansatzYl 



i°° axx 2 f a „(x) + cj§. r^^w + cj§. 

dt Jo dt Jo drfM dt Jo dtfv 



poo r- -, roc 

C A o / dxF A0 (x) 2g A0 x - 2kx 2 - C A0 f A / dxx 2 F A0 (x) 
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h A f°° 

Cai / dx A x A F A1 (x A ) 
Jo 



2 

x [c 



( poo poo \ 

\ c bi J Q dx B x 2 B F m (x B ) + C B o J o dx B x 2 B F m (x B )\ . (A.5) 

We have similar five equations for gene B, which are expressed by the exchange of the 
indexes (A <-> B) for equations (1A.1I) ~ (1A.5[) . We note that there are restrictions for 



the normalization of probability C a o = 1 — C a \. 

When we use the gamma distribution f[2"B"j) for the superposition ansatz, the integral 
factors in equations flA.lj) ~ (IA.5[) are simply replaced by 



dF ni (x) d r°° . ,, (2) 











dxx Tp— = — tjy / dxxF ai (x) = fx { J , (A.6) 
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dxx ^^ - (A 
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dxx 2 ^l=2^(^ + (^) 2 ), (A.9) 

OjJ, 



"at 



where a G {A, B} and i G {0, 1}. 

In order to evaluate the time evolution of the parameters related to gene A 
numerically, we need to calculate dC A \/dt, dC A o/dt, dfj, A l/dt, dfi A l/dt, dfi^l/dt, and 
dfi^l/dt. From equation (IA.1I) . we have dC A i/dt, and then dC AQ /dt is calculated using 

dC A0 dC A \ 

dt ~ d~T' ^ ' ' 

Because equations ( 1A.2I) and (1A.4I) are linear simultaneous equations in dfi A \/dt and 
dfi A \/dt, it is easy to calculate dfi A \/dt and dfi A l/dt. dfi^l/dt and dfi^l/dt are also 
calculated from linear simultaneous equations ( 1A.3I) and ( 1A.5I) . For the time evolution 
of the parameters related to gene B, we perform the same procedures. 
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